function shooting2

%  Solves the BVP using shooting with RK4
%       y'' = f(x,y,y')   for xL < x < xr  '

%  where
%      y(xl) = yL  and y(xR) = yR

% clear all previous variables and plots
clear *
clf

% set boundary conditions and parameters
	xL=0; yL=1;
	xR=1; yR=1;
	
% define title and axes used in plot
hold on
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
title('BVP: Shooting with RK4 and Secant Method','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
axis([0 1 0 1]);
set(gca,'ytick',[0 0.2 0.4 0.6 0.8 1.0]);
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 

% loop used to increase N value
for in=1:3

	% set number of points along the x-axis
	if in==1
		n=10;
	elseif in==2
		n=20;
	elseif in==3
		n=120;
	end

	% generate the points along the x-axis, x(1)=xL and x(n+2)=xR
	x=linspace(xL,xR,n+2);
	h=x(2)-x(1);

	% calculate and then plot solution
	sa=(yR-yL)/(xR-xL);
	y0=[yL sa];
	y=rk4('de_f',x,y0,h,n+2);
	ga=y(n+2)-yR;
	if sa==0
		sb=1;
	else
		sb=0.1*sa;
	end;
	y0=[yL sb];
	y=rk4('de_f',x,y0,h,n+2);
	gb=y(n+2)-yR;

	%  secant method
	tol=0.000001;
	counter=0;
	while abs(sa-sb)>tol
		counter=counter+1;
		sc=sb-(sb-sa)*gb/(gb-ga);
		y0=[yL sc];
		y=rk4('de_f',x,y0,h,n+2);
		gc=y(n+2)-yR;
		sa=sb;  ga=gb;
		sb=sc;  gb=gc;
	end;

	error1=y(n+2)-yR;
	fprintf('\n  n = %3.0f\n  Secant Iterations = %2.0f,  Error at x = 1: %7.2e\n',n,counter,error1) 
	
	% plot the solution
	if in==1
		plot(x,y,'--r','LineWidth',1,'MarkerSize',7)
		legend(' N = 10',3);
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
		pause
	elseif in==2
		plot(x,y,'--b','LineWidth',1,'MarkerSize',7)
		legend(' N = 10',' N = 20',3);
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
		pause
	elseif in==3
		plot(x,y,'--m','LineWidth',1,'MarkerSize',7)
		legend(' N = 10',' N = 20',' N = 120',3);
	end	
	% Set legend font to 14/bold    
	set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
end
hold off

function g=q(x)
ep=0.01; g=-1/ep;

function g=p(x)
ep=0.01; g=-x^2/ep;

function g=f(x)
g=0;


% right-hand side of DE
function z=de_f(x,y)
z=[y(2) -p(x)*y(2)-q(x)*y(1)+f(x)];

% RK4 method
function ypoints=rk4(f,x,y0,h,n)
y=y0;
ypoints=y0(1);
for i=2:n
	k1=h*feval(f,x(i-1),y);
	k2=h*feval(f,x(i-1)+0.5*h,y+0.5*k1);
	k3=h*feval(f,x(i-1)+0.5*h,y+0.5*k2);
	k4=h*feval(f,x(i),y+k3);
	yy=y+(k1+2*k2+2*k3+k4)/6;
	ypoints=[ypoints, yy(1)];
	y=yy;
end;